/**
 * 
 */
package ua.lviv.franko.mlsFunctions.twoDimension;

import java.util.ArrayList;
import java.util.List;

import org.jzy3d.maths.Range;
import org.jzy3d.plot3d.builder.Mapper;

import ua.lviv.franko.integration.IFEMFunction2D;
import ua.lviv.franko.integration.IFunction2D;
import ua.lviv.franko.matrix.SimpleLU;
import ua.lviv.franko.mlsFunctions.oneDimension.weightFunctions.CubicSpline;
import ua.lviv.franko.mlsFunctions.twoDimension.weightFunctions.IWeightFunction2D;
import ua.lviv.franko.mlsFunctions.twoDimension.weightFunctions.SquareSupport;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.fem.LinearRectangularElement;
import ua.lviv.franko.triangulation.SaveStructures.Node;
import ua.lviv.franko.visualisation.Visualise3D;

/**
 * @author marko
 * 
 */
public class BaseFunctionCoupled extends BaseFunction implements IFunction2D {

	protected List<IFEMFunction2D>	femList;

	public BaseFunctionCoupled(int PolinomialDegree, Particle xi, IWeightFunction2D func, List<IFEMFunction2D> femList) {
		super(PolinomialDegree, xi, func);
		this.femList = femList;
	}

	public double calculate(double x, double y) {
		// if((xi.xi - xi.rho>x)||(xi.xi + xi.rho<x)){
		// return 0;
		// }
		// if((xi.yi - xi.rho>y)||(xi.yi + xi.rho<y)){
		// return 0;
		// }
		ArrayList<Double> BaseVector = getBaseVector(x, y, xi);
		int n = BaseVector.size();
		b = new double[n];
		// формуємо матрицю Грамма М
		A = new double[n][n];
		for (int i = 0; i < NodeList.size(); i++) {
			Particle particle = NodeList.get(i);
			calculateAlfa(x, y, particle);
		}

		double val = func.calculate(x, y, xi.xi, xi.yi, xi.rho);
		BaseVector = getBaseVector(xi.xi, xi.yi, xi);
		boolean key = false;
		for (int i = 0; i < b.length; i++) {
			b[i] = BaseVector.get(i) * val;
			if (Math.abs(b[i]) > EPS) {
				key = true;
			}
		}

		double returnValue = 0.0;
		if (key) {
			SimpleLU lu = new SimpleLU(A);
			double[] res = lu.calculate(b);

			BaseVector = getBaseVector(x, y, new Particle(x, y, 0.0));
			double[] femInfluence = this.calculateFEM(x, y);
			for (int i = 0; i < b.length; i++) {
				returnValue += (BaseVector.get(i) - femInfluence[i]) * res[i];
			}
		}
		if (Double.isNaN(returnValue)) {
			return 0;
		}
		return returnValue;
	}

	public double calculateDerivatyX(double x, double y) {
		// if((xi.xi - xi.rho>x)||(xi.xi + xi.rho<x)){
		// return 0;
		// }
		ArrayList<Double> BaseVector = getBaseVector(x, y, xi);
		int n = BaseVector.size();
		b = new double[n];
		// формуємо матрицю Грамма М
		A = new double[n][n];
		for (int i = 0; i < NodeList.size(); i++) {
			Particle particle = NodeList.get(i);
			calculateAlfa(x, y, particle);
		}

		double val = func.calculate(x, y, xi.xi, xi.yi, xi.rho);// func.calculate(x,
																// xi.xi,
																// xi.rho);
		BaseVector = getBaseVector(x, y, new Particle(x, y, 0.0));
		boolean key = false;
		for (int i = 0; i < b.length; i++) {
			b[i] = BaseVector.get(i) * val;
			if (Math.abs(b[i]) > EPS) {
				key = true;
			}
		}

		double returnValue = 0.0;
		if (key) {
			SimpleLU lu = new SimpleLU(A);
			double[] res = lu.calculate(b);

			// val = func.evaluateDerivaty((x-xii)/rho);
			key = false;
			double[] tmp = new double[b.length];
			for (int j = 0; j < this.NodeList.size(); j++) { // по вузлах
				Particle xj = NodeList.get(j);
				val = func.calculateDerivatyX(x, y, xj.xi, xj.yi, xj.rho);// func.calculateDerivaty(x,
																			// xj.xi,
																			// xj.rho);//
																			// func.evaluateDerivaty((x-xj.x)/xj.rho);
				BaseVector = this.getBaseVector(xj.xi, xj.yi, xj);
				double tmpVal = 0.0;
				for (int k = 0; k < res.length; k++) {// альфа на вектор Р і
														// вагу ...
					tmpVal += BaseVector.get(k) * res[k];
				}

				for (int k = 0; k < res.length; k++) {// число на вектор Р
					tmp[k] += BaseVector.get(k) * tmpVal * val;
				}
			}
			val = func.calculateDerivatyX(x, y, xi.xi, xi.yi, xi.rho);// func.calculateDerivaty(x,
																		// xi.xi,
																		// xi.rho);
																		// //func.evaluateDerivaty((x-xii)/rho);
			BaseVector = this.getBaseVector(xi.xi, xi.yi, xi);
			for (int i = 0; i < b.length; i++) {
				b[i] = BaseVector.get(i) * val - tmp[i];
				if (Math.abs(b[i]) > EPS) {
					key = true;
				}
			}
			if (key) {
				double[] res2 = lu.calculate(b);
				BaseVector = this.getBaseVectorDerivatyX(x, y, new Particle(x, y, 0.0));
				for (int j = 0; j < res2.length; j++) {
					returnValue += BaseVector.get(j) * res[j];
				}
				BaseVector = this.getBaseVector(x, y, new Particle(x, y, 0.0));
				for (int j = 0; j < res2.length; j++) {
					returnValue += BaseVector.get(j) * res2[j];
				}
			}

		}
		if (Double.isNaN(returnValue)) {
			return 0;
		}
		// double h=0.01;
		return returnValue;// ((this.calculate(x, y) - this.calculate(x-h, y)
							// )/h);
		// return 0.0;

	}

	public double[] calculateDerivatyXFEM(double x, double y) {
		ArrayList<Double> baseVector = this.getBaseVector(x, y, xi);
		double[] res = new double[baseVector.size()];
		double tmp = 0.0;
		for (int i = 0; i < femList.size(); i++) {
			IFEMFunction2D fem = femList.get(i);
			double xi = fem.getXi();
			double yi = fem.getYi();
			baseVector = this.getBaseVector(xi, yi, new Particle(xi, yi, 0.0));
			tmp = fem.calculateDerivatyX(x, y);
			for (int j = 0; j < res.length; j++) {
				res[j] += baseVector.get(j) * tmp;
			}
		}
		return res;
	}

	public double calculateDerivatyY(double x, double y) {
		// if((xi.xi - xi.rho>x)||(xi.xi + xi.rho<x)){
		// return 0;
		// }
		ArrayList<Double> BaseVector = getBaseVector(x, y, xi);
		int n = BaseVector.size();
		b = new double[n];
		// формуємо матрицю Грамма М
		A = new double[n][n];
		for (int i = 0; i < NodeList.size(); i++) {
			Particle particle = NodeList.get(i);
			calculateAlfa(x, y, particle);
		}

		double val = func.calculate(x, y, xi.xi, xi.yi, xi.rho);// func.calculate(x,
																// xi.xi,
																// xi.rho);
		BaseVector = getBaseVector(x, y, new Particle(x, y, 0.0));
		boolean key = false;
		for (int i = 0; i < b.length; i++) {
			b[i] = BaseVector.get(i) * val;
			if (Math.abs(b[i]) > EPS) {
				key = true;
			}
		}

		double returnValue = 0.0;
		if (key) {
			SimpleLU lu = new SimpleLU(A);
			double[] res = lu.calculate(b);

			// val = func.evaluateDerivaty((x-xii)/rho);
			key = false;
			double[] tmp = new double[b.length];
			for (int j = 0; j < this.NodeList.size(); j++) { // по вузлах
				Particle xj = NodeList.get(j);
				val = func.calculateDerivatyY(x, y, xj.xi, xj.yi, xj.rho);// func.calculateDerivaty(x,
																			// xj.xi,
																			// xj.rho);//
																			// func.evaluateDerivaty((x-xj.x)/xj.rho);
				BaseVector = this.getBaseVector(xj.xi, xj.yi, xj);
				double tmpVal = 0.0;
				for (int k = 0; k < res.length; k++) {// альфа на вектор Р і
														// вагу ...
					tmpVal += BaseVector.get(k) * res[k];
				}

				for (int k = 0; k < res.length; k++) {// число на вектор Р
					tmp[k] += BaseVector.get(k) * tmpVal * val;
				}
			}
			val = func.calculateDerivatyY(x, y, xi.xi, xi.yi, xi.rho);// func.calculateDerivaty(x,
																		// xi.xi,
																		// xi.rho);
																		// //func.evaluateDerivaty((x-xii)/rho);
			BaseVector = this.getBaseVector(xi.xi, xi.yi, xi);
			for (int i = 0; i < b.length; i++) {
				b[i] = BaseVector.get(i) * val - tmp[i];
				if (Math.abs(b[i]) > EPS) {
					key = true;
				}
			}
			if (key) {
				double[] res2 = lu.calculate(b);
				BaseVector = this.getBaseVectorDerivatyY(x, y, new Particle(x, y, 0.0));
				for (int j = 0; j < res2.length; j++) {
					returnValue += BaseVector.get(j) * res[j];
				}
				BaseVector = this.getBaseVector(x, y, new Particle(x, y, 0.0));
				for (int j = 0; j < res2.length; j++) {
					returnValue += BaseVector.get(j) * res2[j];
				}
			}

		}
		if (Double.isNaN(returnValue)) {
			return 0;
		}
		return returnValue;
		// return 0.0;

	}

	public double[] calculateDerivatyYFEM(double x, double y) {
		ArrayList<Double> baseVector = this.getBaseVector(x, y, xi);
		double[] res = new double[baseVector.size()];
		double tmp = 0.0;
		for (int i = 0; i < femList.size(); i++) {
			IFEMFunction2D fem = femList.get(i);
			double xi = fem.getXi();
			double yi = fem.getYi();
			baseVector = this.getBaseVector(xi, yi, new Particle(xi, yi, 0.0));
			tmp = fem.calculateDerivatyY(x, y);
			for (int j = 0; j < res.length; j++) {
				res[j] += baseVector.get(j) * tmp;
			}
		}
		return res;
	}

	public double[] calculateFEM(double x, double y) {
		ArrayList<Double> baseVector = this.getBaseVector(x, y, xi);
		double[] res = new double[baseVector.size()];
		double tmp = 0.0;
		for (int i = 0; i < femList.size(); i++) {
			IFEMFunction2D fem = femList.get(i);
			double xi = fem.getXi();
			double yi = fem.getYi();
			baseVector = this.getBaseVector(xi, yi, new Particle(xi, yi, 0.0));
			tmp = fem.calculate(x, y);
			for (int j = 0; j < res.length; j++) {
				res[j] += baseVector.get(j) * tmp;
			}
		}
		return res;
	}

	/**
	 * @param args
	 */
	public static void main(String[] args) {
		// final LinearRectangularElement element = new LinearRectangularElement(new Node(0, 0, true, 0), new Node(0.5, 0.0, true, 1), new Node(0.5, 0.5, true, 2), new Node(0, 0.5,
		final LinearRectangularElement element = new LinearRectangularElement(new Node(0, 0, true, 0), new Node(1.0, 0.0, true, 1), new Node(1.0, 1.0, true, 2), new Node(0, 1.0,
				true, 3), 1, 1, 1, null, null);
		List<IFEMFunction2D> fem = new ArrayList<IFEMFunction2D>();
		fem.add(new IFEMFunction2D() {

			@Override
			public double calculateDerivatyY(double x, double y) {
				return 0;
			}

			@Override
			public double calculateDerivatyX(double x, double y) {
				return 0;
			}

			@Override
			public double calculate(double x, double y) {
				return element.N1.calculate(x, y);
			}

			@Override
			public double getYi() {
				return element.getA().y;
			}

			@Override
			public double getXi() {
				return element.getA().x;
			}
			@Override
			public Node getNode() {
				return element.getA();
			}
		});
		fem.add(new IFEMFunction2D() {

			@Override
			public double calculateDerivatyY(double x, double y) {
				return 0;
			}

			@Override
			public double calculateDerivatyX(double x, double y) {
				return 0;
			}

			@Override
			public double calculate(double x, double y) {
				return element.N2.calculate(x, y);
			}

			@Override
			public double getYi() {
				return element.getB().y;
			}

			@Override
			public double getXi() {
				return element.getB().x;
			}
			@Override
			public Node getNode() {
				return element.getB();
			}
		});
//		fem.add(new IFEMFunction2D() {
//
//			@Override
//			public double calculateDerivatyY(double x, double y) {
//				return 0;
//			}
//
//			@Override
//			public double calculateDerivatyX(double x, double y) {
//				return 0;
//			}
//
//			@Override
//			public double calculate(double x, double y) {
//				return element.N3.calculate(x, y);
//			}
//
//			@Override
//			public double getYi() {
//				return element.getC().y;
//			}
//
//			@Override
//			public double getXi() {
//				return element.getC().x;
//			}
//		});
		fem.add(new IFEMFunction2D() {

			@Override
			public double calculateDerivatyY(double x, double y) {
				return 0;
			}

			@Override
			public double calculateDerivatyX(double x, double y) {
				return 0;
			}

			@Override
			public double calculate(double x, double y) {
				return element.N4.calculate(x, y);
			}

			@Override
			public double getYi() {
				return element.getD().y;
			}

			@Override
			public double getXi() {
				return element.getD().x;
			}
			@Override
			public Node getNode() {
				return element.getD();
			}
		});
		ArrayList<Particle> nodes = new ArrayList<Particle>();
		// параметри області:
		double ax = 0, bx = 1, ay = 0, by = 1;

		int n = 10, m = 10;

		double hx = (bx - ax) / (double) n, hy = (by - ay) / (double) m;

		// вузли+функції:
		double rho = 2.5 * hx;
		for (int i = 1; i < n; i++) {
			for (int j = 1; j < m; j++) {
				Particle p = new Particle(ax + i * hx, ay + j * hy, rho);
				nodes.add(p);
			}
		}
		final BaseFunctionCoupled ff = new BaseFunctionCoupled(1, nodes.get(0), new SquareSupport(new CubicSpline()), fem);
		ff.setNodeList(nodes);
//		ff.calculate(0.5, 0.5);
//		ff.calculate(0.75, 0.75);

		// Mapper mapper2 = new Mapper() {
		//
		// @Override
		// public double f(double x, double y) {
		// double h = 0.01;
		// double tmp = (ff.calculate(x, y) - ff.calculate(x-h, y) )/h;
		// return tmp;//(x, y);
		// }
		// };

		Visualise3D.show(new Mapper() {
			@Override
			public double f(double x, double y) {
				// return ff.calculateDerivatyX(x, y);//(x, y);
				// double val =ff.calculateDerivatyX(x, y);
				// double h = 0.00001;
				long start = System.currentTimeMillis();
				double val = ff.calculate(x, y);
				System.err.println(System.currentTimeMillis()-start);
				// double val1 = (val - ff.calculate(x - h, y)) / h;
				// double val2 = (val - ff.calculate(x, y - h)) / h;
				// System.out.println(x + "\t" + y + "\t" + val + "\t" + val1 + "\t" + val2);
				return val;
			}
		}, new Range(0.0, 1.0), new Range(0.0, 1.0));

	}

}
